#from numpy import *
from mpl_toolkits.mplot3d import Axes3D
from matplotlib import cm
from pylab import * # needed for 2D plot
import matplotlib.pyplot as plt
import numpy as np

data_L1 = np.genfromtxt("out_L1.txt");
data_L2 = np.genfromtxt("out_L2.txt");

t_L1 = data_L1[:,0] 
e_L1 = data_L1[:,1]
r_L1 = data_L1[:,2]

t_L2 = data_L2[:,0] 
e_L2 = data_L2[:,1]
r_L2 = data_L2[:,2]

theta_L1 = np.radians(t_L1)
ele_L1   = np.radians(e_L1)
h_L1     = r_L1 * np.cos(ele_L1) # work out the horizontal component first

theta_L2 = np.radians(t_L2)
ele_L2   = np.radians(e_L2)
h_L2     = r_L2 * np.cos(ele_L2) # work out the horizontal component first

N_L1 = (h_L1) * np.cos(theta_L1)
E_L1 = (h_L1) * np.sin(theta_L1)
U_L1 = (r_L1) * np.sin(ele_L1)

N_L2 = (h_L2) * np.cos(theta_L2)
E_L2 = (h_L2) * np.sin(theta_L2)
U_L2 = (r_L2) * np.sin(ele_L2)

#=============================================================
# 3D plot of PCV corrections in NEU 
#=============================================================
fig = plt.figure(1)
ax = Axes3D(fig)
ax.scatter(N_L1, E_L1, U_L1, N_L1, color='r')
#ax.scatter(N_L2, E_L2, U_L2, N_L2, color='b')
title('PCV coreection vs Elevation and Azimuth')

#=============================================================
# Plot all of the L1 & L2 PCV corrections per elevation angle 
#=============================================================
##plt.figure(2)
#plt.subplot(211)
# Try just a 2d plot
##r_L1.shape = 72,18
#r_L2.shape = 72,19
#m = linspace(0,90,19)
##m = linspace(5,90,18)
##(l,j) = r_L1.shape
##for i in range(0, l):
##    plt.plot(m, r_L1[i,:], 'r.' )
#    plt.plot(m, r_L2[i,:], 'b.' )

##print m
##title('PCV correction vs Elevation')
##xlabel('Elevation in (degrees)')
##ylabel('PCV correction (mm)')
#=============================================================
# Plot the mean PCV correction versus elevation angle
#=============================================================
##plt.figure(3)
#for i in range(0,19):
##for i in range(1,18):
##    m_L1 = r_L1[:,i].mean()
##    s_L1 = r_L1[:,i].std() * 3
#    m_L2 = r_L2[:,i].mean()
#    s_L2 = r_L2[:,i].std() * 3
##    e = i * 5
##    errorbar(e,m_L1,s_L1,fmt='ro')
#    errorbar(e,m_L2,s_L2,fmt='bo')

##grid(True)
##title('Mean PCV correction vs Elevation')
##xlabel('Elevation (degrees)')
##ylabel('PCV correction (mm)')

#=============================================================
# 3D plot of PCV corrections in AZ,El,R 
#=============================================================
##fig = plt.figure(4)
##ax = Axes3D(fig)
##z_L1 = data_L1[:,2]
#z_L2 = data_L2[:,2]
##ax.scatter(t_L1, e_L1, z_L1, t_L1, color='r')
#ax.scatter(t_L2, e_L2, z_L2, t_L2, color='b')

##fig = plt.figure(5)
##ax = Axes3D(fig)

##X = outer( ones(size(e_L1)), t_L1 ) #+ 0.37
##Y = outer( ones(size(e_L1)), e_L1 ) # + 0.86
##Z = outer( ones(size(e_L1)), data_L1[:,2]  )# + 90.02
##ax.plot_wireframe(X, Y, Z, rstride=1, cstride=1, cmap=cm.jet)
#ax.w_zaxis.set_major_locator(LinearLocator(6))


#fig = plt.figure(6)
#ax = Axes3D(fig)

#X = outer( ones(size(e_L1)), t_L1 ) #+ 0.37
#Y = outer( ones(size(e_L1)), e_L1 ) # + 0.86
#Z = outer( ones(size(e_L1)), data_L1[:,2]  )# + 90.02
#ax.plot_surface(X, Y, Z, rstride=5, cstride=5, cmap=cm.jet, linewidth=5, antialiased=False)
plt.show()



